PAM-664 


THE  construction  OF  ORTHOGONAL  EIGENVECTORS 
FOR  TIGHT  CLUSTERS  BY  USE  OF  SUBMATRICES 


Beresford  Parlett 


Department  of  Mathematics 
University  of  California 
Berkeley,  California  94720 


ii' ;  : Uainputadv-:  _ 


1997071/  157 


PTIC  QUALET^  iHSFSCEl 


This  report  was  done  with  support  from  the  Center  for  Pure  and 
Applied  Mathematics.  Any  conclusions  or  opinions  in  this  report  repre¬ 
sent  solely  those  of  the  author(s)  and  not  necessarily  those  of  the  Center 
for  Pure  and  Applied  Mathematics  or  the  Department  of  Mathematics. 


DEPARTMENT  OF  THE  NAVY 

OFFICE  OF  NAVAL  RESEARCH 
SEATTLE  REGIONAL  OFFICE 
1107  NE  45TH  STREET.  SUITE  350 

SEATTLE  WA  98105-4631  in  reply  refer  tq; 


4330 
ONR  247 
1 1  Jul  97 


From:  Director,  Office  of  Naval  Research,  Seattle  Regional  Office,  1 107  NE  45th  St.,  Suite  350, 
Seattle,  WA  98105 

To:  Defense  Technical  Center,  Attn:  P.  Mawby,  8725  John  J.  Kingman  Rd.,  Suite  0944, 

Ft.  Belvoir,  V A  22060-62 18 

Subj:  RETURNED  GRANTEE/CONTRACTOR  TECHNICAL  REPORTS 


1.  This  confirms  our  conversations  of  27  Feb  97  and  1 1  Jul  97.  Enclosed  are  a  number  of 
technical  reports  which  were  returned  to  our  agency  for  lack  of  clear  distribution  availability 
statement.  This  confirms  that  all  reports  are  unclassified  and  are  “APPROVED  FOR  PUBLIC 
RELEASE”  with  no  restrictions. 


2.  Please  contact  me  if  you  require  additional  information.  My  e-mail  is  silverr@,onr.  navy,  mil 
and  my  phone  is  (206)  625-3196. 


/ 


/ 


ROBERT  J.  SILVERMAN 


The  Construction  of  Orthogonal  Eigenvectors 
for  Tight  Clusters  by  Use  of  Submatrices 

Beresford  N.  Parlett  * 

January  18,  1996 


Abstract 

The  goal  is  to  compute  eigenvectors  of  a  symmetric  tridiagonal 
matrix  T  that  are  orthogonal  to  working  accuracy.  Consider  a  cluster 
of  Tti  very  close  eigenvalues  that  are  reasonably  well  separated  from 
the  remaining  spectrum.  We  show  here  that  there  are  m  principal 
submatrices  of  T  such  that  only  the  nearest  neighbors  overlap  and 
each  submatrix  has  a  simple,  isolated  eigenvalue  in  the  convex  huU 
of  the  cluster  with  eigenvector  having  small  entries  in  the  first  and 
last  positions.  This  eigenvector  is  padded  with  zero  entries,  above 
and  below,  to  make  it  conform  to  T.  The  set  of  vectors,  one  from 
each  submatrix,  forms  a  good  basis  for  the  invariant  subspace.  Each 
basis  vector  may  be  modified,  if  necessary,  by  its  nearest  neighbors  to 
produce  an  orthonormal  basis. 

The  only  communication  that  may  be  needed,  in  such  situations, 
is  between  nearest  neighbors. 

We  give  a  good  bound  on  the  dot  product  of  nearest  neighbors.  A 
variety  of  examples  illustrate  the  theory. 

The  ideas  in  this  paper  were  presented  at  the  ENUMATH  meeting  at 
CNRS,  Paris,  in  September  1995  and  at  the  ILAY  workshop  at  Cerfacs  in 
October,  1995. 
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versity  of  California,  Berkeley,  CA  94720,  USA.  Supported  by  ONR,  Contract  N000014- 
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1  Introduction 

^Inverse  Iteration  gives  a  very  satisfactory  solution  to  the 
problem  as  far  as  reasonably  well  separated  eigenvalues  are  con¬ 
cerned.  The  problem  of  determining  reliably  full  digital  infor¬ 
mation  in  the  subspace  spanned  by  eigenvectors  corresponding 
to  coincident  or  pathologically  close  eigenvalues  has  never  been 
satisfactorily  solved.’ 

J.  H.  Wilkinson 

(from  ‘The  Algebraic  Eigenvalue  Problem’,  1965, 

Chapter  5,  p.  344.) 

This  quotation  is  over  30  years  old  and  yet  most  experts  would  agree  that 
its  claim  is  still  true  in  1995. 

Our  goal  is  to  build  up  an  orthonormal  basis  for  the  invariant  subspace 
associated  with  a  reasonably  isolated  cluster  of  very  close  eigenvalues  of  a 
symmetric  tridiagonal  matrix  T.  This  is  achieved  b}'  the  QR  algorithm  and 
it  is  only  the  relatively  high  cost  of  accumulating  all  the  plane  rotations  that 
drives  the  search  for  other  techniques.  Certain  codes  (e.  g.,  xstein)  currently 
used  in  LAPACK,  and  other  libraries  such  as  NAG  and  IMSL  and  ESSL, 
decline  in  both  efficiency  and  quality  of  output  as  eigenvalues  get  closer  to 
each  other.  The  reason  is  that  these  codes  are  based  on  inverse  iteration  and 
it  is  extremely  difficult  to  choose  automatically  suitable  right  hand  sides  to 
ensure  both  the  spanning  property  (accuracy)  and  orthogonality. 

One  way  out  of  the  difficulty  is  to  discard  an  appropriate  set  of  rows  from 
the  top  and  bottom  of  T  and  to  work  with  the  remaining  submatrix  to  obtain 
a  basis  vector.  The  well  known  test  matrix  VV^i  that  is  discussed  in  Section  5 
has  its  largest  two  eigenvalues  equal  to  single  (and  double)  precision.  In  this 
easy  case  it  suffices  to  compute  the  eigenvectors  of  the  submatrices  in  rows 
1  to  19  and  3  to  21,  append  zero  entries  at  the  bottom  of  one  and  the  top  of 
the  other,  and  finally  deliver  the  internal  and  external  bisectors  of  these  two 
vectors.  Even  the  tiny  entries  are  computed  to  high  relative  accuracy. 

The  guiding  principle  behind  this  approach  is  that  an  eigenvector  associ¬ 
ated  with  a  simple,  isolated  eigenvalue  is  easy  to  compute.  So,  for  a  cluster 
of  m  close  eigenvalues,  the  task  is  to  find  m  different  submatrices  of  T  each 
of  which  has  an  isolated  eigenvalue  in  the  convex  hull  of  the  cluster.  That  is 
not  enough.  The  eigenvector  of  any  internal  submatrix  must  have  small  en- 
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tries  in  its  first  and  last  positions.  It  is  not  obvious  that  these  specifications 
can  always  be  met  and  the  goal  of  this  paper  is  to  provide  the  theory  that 
supports  our  use  of  submatrices. 

It  turns  out  that  the  case  of  close  pairs  is  fundamental.  The  general  case 
may  be  reduced  to  considering  a  number  of  pairs.  A  surprising  outcome  of 
these  investigations  is  that  the  difficulty  of  computing  an  accurate  orthogonal 
basis  is  not  properly  measured  by  the  gap  between  eigenvalues.  The  support 
of  the  vectors  (the  positions  holding  nonnegligible  values)  plays  an  important 
role.  Strictly  speaking  zero  entries  in  eigenvectors  are  extremely  rare  and, 
when  they  do  occur,  they  are  isolated.  Since  we  ignore  isolated  zero  entries 
we  might  conclude  that  the  support  of  all  the  eigenvectors  we  seek  is  the  full 
index  set.  To  take  such  a  pedantic  view  would  be  to  miss  an  important  fact 
about  many,  but  not  all,  n  x  n  tridiagonal  matrices  as  n  grows:  the  active 
part  of  the  eigenvector  is  confined  to  a  small  part  of  the  domain,  perhaps 
only  30  or  40  consecutive  positions.  The  remaining  entries  carry  a  tiny  bit 
of  noise.  Consequently  we  use  the  term  support  somewhat  informally.  If  the 
support  of  two  eigenvectors  is  disjoint  then  they  are  orthogonal  however  close 
the  associated  eigenvalues  may  be. 

Closely  related  to  the  idea  of  support  is  the  concept  of  the  overlap  of  two 
vectors  x  and  y; 

OvcTlap(x,y)  :=  |j^]|T|j^. 

the  cosine  between  the  vectors  of  absolute  values.  Our  two  main  results  are: 
1.  Each  normalized  vector  x  produced  from  an  appropriate  submatrix  by 
appending  zeros  has  p  :=  x*Tx  in  the  cluster  interval  and 


||(r  -  pl)x\\  =  0{clusitr  length) 

2.  Any  two  normalized  vectors  x  and  y  from  overlapping  submatrices  satisfy 

// cluster  length\^^‘^ 

Overlap{x,y)  =  0  ( - —  j 


VV 


gap 


where  gap  is  the  separation  of  the  cluster  from  the  remaining  spectra  of  the 
two  submatrices. 

Submatrices  only  overlap  their  nearest  neighbors.  Physicists  would  say 
that  the  overlap  matrix  is  tridiagonal  and  close  to  the  identity.  However 
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the  subspace  spanned  by  these  subvectors  is  more  accurate  than  indicated 
by  the  residual  norms  of  the  subvectors  separately.  That  is  the  content  of 
Section  9.  What  Results  1  and  2  show  is  that  the  subvectors  yield  a  sparse, 
nearly  orthogonal  basis  for  a  good  subspace. 

The  results  developed  here  are  pure  matrix  theory,  there  is  no  reference 
to  machine  precision.  They  depend  strongly  on  the  tridiagonal  from. 

The  first  result  is  somewhat  surprising  because  experience  with  the  Lan- 
czos  algorithm  had  suggested  that  the  best  that  could  be  guaranteed  for  any 
vector  X  obtained  from  a  submatrix  was  \\{T  —  fil)x\\  =  0{\/ cluster  length) 
and  that  is  not  strong  enough  for  our  purposes. 

Tight  clusters  of  close  eigenvalues  are  not  the  only  spectral  distributions 
to  pose  challenges  to  numerical  analysts.  Consider  eigenvalues  that  form  a 
geometric  progression.  This  case  is  not  troublesome  if  the  eigenvalues  are 
computed  to  high  relative  relative  accuracy.  A  different  problem  is  posed 
by  small  perturbations  of  the  identity  matrix  I.  How  large  must  the  per¬ 
turbation  be  before  the  user  will  not  accept  the  columns  of  I  as  suitable 
eigenvectors?  Shifting  by  I  and  scaling  by  the  largest  remaining  entry  makes 
the  perturbations  seem  important  and  the  shifted  and  scaled  eigenvalues 
seem  well  separated.  One  would  then  compute  a  set  of  eigenvectors  very 
different  from  the  columns  of  7.  It  would  seem  that  only  the  user  can  choose 
between  the  two  and  we  shall  not  pursue  this  question  here. 

We  begin  our  analysis  in  Section  6  by  introducing  the  envelope  of  a  cluster 
and  showing  how  it  guarantees  Result  1  with  a  constant  bounded  by  ^njl 
but  normally  much  smaller.  The  envelope  also  reveals  good  choices  for  sub¬ 
matrices.  Then  we  take  a  different  approach  to  show  how  the  submatrix 
indices  affect  the  constants  behind  the  0{cluster  length).  Our  results  make 
heavy  use  of  detailed  properties  of  tridiagonal  matrices  and  that  material  is 
gathered  in  Section  3  and  Section  10. 

The  reader  is  urged  to  read  Section  2  on  notation. 


2  Notation 


Pi  P2  •  Pn-l 

ai  02  •  •  ct„ 

Pi  P2  •  Pn-1 


T  =  tridiag 


P2 
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Eigenvalues:  Ai  <  A2  <  •  •  •  ^  A„.  (1) 

Normalized  eigenvectors:  81,82,  ,Sn-  (2) 

When  /3,-  ^  0,  i  =  1, . . . ,  n  -  1  then  T  is  unreduced.  In  that  case  the  inequal¬ 
ities  in  (1)  are  strict. 

The  principal  submatrix  of  T  in  rows  j,j  -H  1»  •  •  •  >  ^  is  denoted  by  either 
'j'j'k  Qj.^  simply,  by  (j  :  k).  The  eigenvalues  of  a  submatrix  (j  :  fc)  are  called 

Ritz  values  or  R- values: 


<0^=''<...<CVr  (3) 

The  normalized  eigenvector  of  is  The  mth  entry  is  sf{m).  The 
natural  way  to  index  the  entries  of  sf  is  from  j  to  k,  not  from  1  to  k  —  j  +  1. 
In  this  way  sf'*'  is  embedded  in  a  vector  conformable  to  T  and  we  consider 
entries  in  positions  1  :  j  —  1  and  A:  -|- 1  :  n  as  zero.  Sometimes  we  write  0^  as 
To  simplify  further  we  write  6^^^  for  when  ^  lies  in  [A_,A.j.],  our 
cluster. 

The  characteristic  polynomial  of  T^'^,  or  {j  :  k),  is  defined  by 

:=  det  {rl  — 

We  write  Xi  for  -  fri  general  column  vectors  are  denoted  by  lower  case  Ro¬ 
man  letters  in  boldface  type:  v,z,. . .,  and  individual  entries  are  v{j),  z{l),  — 
The  size  of  the  identity  matrix  I  is  given  by  the  context  and  its  columns  are 

e,-  =  (0, ...  0, 1, 0  ... ,  0)^  1  in  position  i. 

3  Eigenvectors  and  Error  Estimates 

The  unreduced  n  x  n  matrix  T  is  uniquely  determined  (up  to  signs)  by 
its  spectrum  {Aa,...,A„}  and  (n  -  1)  appropriate  extra  items.  This  extra 
information  may  take  various  forms:  the  spectrum  of  submatrix  (1  :  n  —  1) 
or  {2  in),  or  the  squares  of  the  top  (or  bottom)  entries  of  the  normalized 
eigenvectors.  See  [4]  or  [2]. 

Consequently  there  are  numerous  expressions  for  the  eigenvectors  of  T 
and  we  give  some  of  them  here. 

Here  we  assume  that  T  is  unreduced:  /?,■  7^  0,  f  =  1, . . . ,  n  -  1- 
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Perhaps  the  oldest  formulae  for  the  unnormalized  eigenvector  for  A  are 

V  ’  /3l  ’  Plh  ’ 

(  X^="(A)  X^-(A) 

/3n-l  ’7 

The  catch  to  using  these  formulae  is  that  they  can  overflow  very  easily  and 
they  may  sometimes  be  extremely  sensitive  to  tiny  changes  in  A. 

There  are  attractive  formulae  for  the  magnitudes  of  the  entities  of  nor¬ 
malized  eigenvectors  Si, . . . ,  Sn-  For  each  i  =  1, 2, . . . ,  n,  and  x  —  X^’” 


s,(l)s,(n)x'(A,)  =  (6) 

^.•(I)^X'(A.)  =  X^=”(A.),  (7) 

5.(n)"x'(A.)  =  X'=”-'(A.),  (8) 

^.■Orx'(A.)  =  x^'-^A.Ox^-'^'-lAO-  (9) 


All  these  formulae  come  from  a  result  in  [4]  that 

adi(A./-r)  =  s,s’x'(A,) 

where  adj{M)  is  the  classical  adjugate  of  M. 

These  results  raise  the  question  of  how  to  find  the  correct  signs  of  the 
entries  and  that  leads  to  the  next  few  observations. 

Given  any  unreduced  t  there  is  a  unique  A  =  diag{±l,  ±1, . . . ,  ±1)  such 
that  T  =  ATA  has  positive  off-diagonal  entries.  Observe  that 

Ts  =  sA  -4=^  r(As)  =  (As)A. 

Consequently  it  is  possible  to  normalize  any  symmetric  T  so  that  it  becomes 
a  direct  sum  of  unreduced  tridiagonals  each  of  which  has  positive  off-diagonal 
entries.  From  now  on  we  assume  that  ^,  >  0,  i  =  1, . . . ,  n  —  1. 

These  Ts  are  oscillation  matrices,  a  term  coined  by  Krein  and  Gantmacher 
[1].  The  number  of  sign  changes  among  consecutive  entries  of  the  eigenvector 
Si  is  n  —  i.  Recall  that  we  label  eigenvalues  so  that 


Ai  <C  A2  <t  .  .  .  <^  An 
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For  example,  for  the  second  difference  SINE  matrix  (1,  2, 1),  the  eigenvalues 
lie  in  -4  <  A,-  <  0  and  the  fundamental  mode  is  given  by  associated  with 
the  rightmost  eigenvalue.  We  use  right  and  left  rather  than  large  or  small  in 
order  to  be  translation  invariant.  For  some  people  it  is  unnatural  to  say  that 
-1  is  larger  than  -2  because  they  think  of  magnitudes  when  the  word  large 
is  used. 

The  correct  signs  may  be  attached  to  magnitudes  when  certain  informa¬ 
tion  is  available.  A  warning  is  in  order  here.  Although  the  characteristic 
polynomial  is  defined  by  x{t)  =  det  (r/  -  T)  most  people  compute  with 
T  -  tI.  For  example,  the  pivots  in  Gaussian  elimination  are  usually  com¬ 
puted  from  the  recurrence 

di+i  =  (a,+i  -  t)  -  ^J/di 

and  so 

/=J(r)  =  {-iydi...dj. 

If  r  =  A,*  (4)  may  be  rewritten  as 


However,  in  practice  r  ^  Xi. 

We  turn  now  to  approximations  and  how  to  assess  them.  Let 

^  l|5fl|  =  l, 

define  an  eigenvector  of  a  submatrix.  Consider  sj  as  an  approximate  eigen¬ 
vector  of  T  and  append  zeros  to  sf*  to  make  it  conform.  Now  drop  the 
subscript  /  and  observe  that 

ll(r  -  +  (11) 

Paige’s  Persistence  Theorem.  Write  for  Olf  For  any  j,  1  <j  <n, 
and  any  i,  1  <i  <  j,  the  closed  interval 

-f  M 
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contains  at  least  one  Ritz  value  fov  each  I  —  1,2,  ...,h  j-  Heve 

Proof.  By  partitioning  T  one  sees  that 

^  ^  =  e,+i/9, 5^0). 

So,  by  a  standard  theorem,  see  Chap.  4  in  [4],  there  is  a  such  that 

|^(i+0  _  <  II  (T^--i+‘  -  )  II  = 

□ 

Frequently  the  bound  ^jUi  is  a  severe  overestimate.  Whenever  is 
isolated  from  the  other  Ritz  values  then  one  defines 

gap{v,j,l)  =  gap{i)  :=  min|0['^'^'^  - 

Gap  Theorem.  For  /  =  1, 2, . . . ,  n  -  j  there  is  a  such  that 
1^0+')  _  6ib)|  <  (fjUPiflgapirJJ). 

See  [4,  Chap.  11]. 

This  is  a  huge  improvement  over  /SjW,-  in  most  cases. 


The  Average-^  Result 


Lemma  1  For  i  =  1, 2, . . . ,  n, 


Proof.  For  each  i,  1  <  i  <  n 


s,(l)  s,(n)  x'(A,)  =  ^1  • . . .  •  /?n-i 


and 
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Since  the  geometric  mean  is  majorized  by  the  arithmetic  mean 
|s,(l)s,(n)|  <  ^  +  <  2' 

Hence  , 

and 

(A  •  ...-A-i )''<”■'*  < 

as  claimed.  P 

Quantities  of  the  form  I3jlgap  will  be  0(1)  on  average  and  they  occur  in 
our  analysis.  Speaking  loosely  we  may  say  that  the  average  jS  is  bounded  by 
the  average  distance  of  any  eigenvalue  from  all  the  others. 


4  On  First  and  Last  Entries 

Let  A_  and  A+  be  two  adjacent  eigenvalues  of  T  well  separated  from  the 
remaining  spectrum.  Let 

rs±  =  S±A±,  ||s±||  =  1. 

From  (6)  in  Section  3, 

5±(l)s±(n)x'(A±)  =  /5i  •  •  •  /5n-l 

Thus 

s_(l)s_(n)  X'('^+) 

(A- — A+)  nj('^-  “  ( TT  \  \  '\ 

=  -nC-^ 


// 

=  -1  +  (A+-A_)  ^  \ 


-1  +  0 


gap  ) 


9 


where  gap  =  minj  |A+  —  Aj|.  So,  when  A4.  —  A_  -C  gap, 


5+(l) 

s_(n) 

s_(l) 

<Nii/ 

s+(n) 

If  |5+(1)|  >  |s_(l)|  then  the  support  of  s+  is  concentrated  in  the  top  of  an 
n- vector  whereas  the  support  of  S-  is  concentrated  in  the  lower  section.  Thus 
the  supports  may  be  nearly  disjoint,  however  close  A_  and  A4.  may  be.  When 
\s+{l)\  w  |s-(l)|  then  s+  and  s_  have  the  same  support  and  it  is  easier  to 
approximate  the  bisectors  {s+  ±  s_)/\/2. 


5  An  Example:  W21 

/  1  1  .  1  \ 

W  =  =  tridiag  I  10  9  •  #9  10  I  . 

V  1  1  •  1  / 

This  matrix,  and  its  companion  W2\,  were  designed  by  Wilkinson,  see  Ch.  6 
in  [6],  to  illustrate  some  subtle  points  concerning  the  computation  of  eigen¬ 
vectors.  Although  the  eigenvalues  of  W  are  distinct  in  exact  arithmetic  there 
are  several  pairs  that  are  equal  to  single  precision  (e  =  1.2  x  10"’’).  However 
the  smaller  eigenvalue  pairs  are  not  so  close  and  the  only  negative  eigenvalue 
is  well  separated  from  the  rest.  Table  1  gives  selected  eigenvalues.  Note  that 
the  W21  is  persymmetric:  invariant  under  reversal. 

Tf  we  separate  W21  into  a  direct  sum  by  putting  (3u  =  0  then  we 
obtain  independent  orthogonal  vectors  which  span  the  subspace 
corresponding  to  the  pathologically  close  pair  of  eigenvalues  A20 
and  A21  to  a  very  high  accuracy.  It  is  therefore  quite  possible  that 
such  a  decomposition  is  always  permissible,  and  what  is  needed  is 
some  reliable  method  of  deciding  when  and  where  to  decompose.’ 

J.  H.  Wilkinson 

(from  ‘The  Algebraic  Eigenvalue  Problem’,  1965, 

Chapter  5,  p.  330.) 

This  quotation  shows  that  Wilkinson  considered  the  possibility  of  using 
disjoint  submatrices  to  obtain  orthogonal  basis  vectors.  Our  investigations 
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A21 

10.746194182903393.. 

A20 

10.746194182903322.. 

Ai3 

6.0002340.. 

Ai2 

6.0002175.. 

A3 

0.9475.. 

A2 

0.2538.. 

Ai 

-1.1254.. 

Table  1:  Selected  Eigenvalues  of 


show  that  this  idea  only  works  in  extreme,  and  easy,  cases.  What  we  illustrate 
here  is  that  the  submatrices  must  be  allowed  to  overlap. 

In  a  nice  recent  study,  see  [7],  Ye  has  turned  Wilkinson  s  idea  into  some 
theorems  relating  a  pair  of  eigenvalues  of  T  to  an  eigenvalue  of  a  certain 
submatrix  and  an  eigenvalue  of  the  complementary  submatrix.  However  he 
does  not  give  results  on  the  related  eigenvectors  and  that  is  our  concern. 

The  largest  pair:  A20  ^^nd  A21  (equal  to  single  precision  at  value  A) 

We  drop  the  last  2  rows  and  approximate  the  Ritz  vector  z+  by  solving 

(H/i:i9_A7)z+ =  ei7i,  1|2+||  =  1- 

We  drop  the  first  2  rows  and  approximate  the  Ritz  vector  z_  by  solving 

=  621721,  Ik-11  =  1- 

We  insert  zero  entries  to  make  z_  and  z+  conform  to  W .  It  turns  out  that 

.  z_  =  1.22  X  10-'^  !,  71  =  -8.33  x  10-^  721  =  -8.33  x  10“^ 

Thus  z_,  and  z+  pass  both  requirements  for  good  eigenvectors,  orthogonal¬ 
ity  and  small  residuals.  However  an  equally  valid  basis  for  the  dominant 
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invariant  2-space  is  {{z+  —  2_)/V2j  (2+  -h  z^)/y/2],  the  bisectors.  It  is 
gratifying  that  this  computed  basis  delivers  the  exact  eigenvectors  correct  to 
single  precision,  even  the  smallest  entries. 

The  results  from  dropping  three  rows  instead  of  two  were  indistinguish¬ 
able  from  the  ones  above.  In  fact  this  case  is  so  easy  that  one  can  use 
submatrices  (1  :  j)  and  (22  —  j  :  21)  for  j  =  14, ...,20,  for  satisfactory 
results. 

The  pair  near  6;  Ai2  and  A13 

The  shift  invariant  measure  of  a  symmetric  matrix  is  its  spread  defined  as 
^max  ~  Amin-  1^1  this  Case 

$pread{W)  =  A21  —  Aj  =  11.87, 
interval  width  =  A13  —  A12  =  11.66  •  c  ■  spread. 

The  two  best  choices  are  (1  :  15),  (7  ;  21)  and  (1  :  17),  (5  :  21).  The  outputs 
from  the  two  choices  are  are  barely  distinguishable  in  single  precision.  The 
Ritz  value  p  for  (1  :  17)  and  (5  :  21)  is  not  exactly  at  the  mean  but,  in  single 
precision,  appears  to  be  so: 


p  =  6.000226. 

Approximate  the  corresponding  Ritz  vectors  by  solving 

=  6676,  -  pI)Z-  =  615715 

in  single  precision  to  find 

-  pl)z.\\  =  1.4  X  10"®  =  10  •  e  •  spread 


and 

•  z_  =  2.1  X  10"®  =  8.4  •  n  •  e. 

Section  9  shows  how  the  subspace  span  {2:-,2+}  is  slightly  more  accurate 
than  indicated  by  the  residual  norms  for  2_  and  z+  separately. 

In  several  ways  {z.,  z+}  is  a  satisfactory  basis  but  the  pair  of  bisectors 
{(z+  —  z_)/\/2,  (z+  -t-  z_)/\/2}  is  even  better.  Just  as  for  the  pair  A20,  A21 
the  computed  bisectors  deliver  the  exact  eigenvectors  to  single  precision,  even 
the  smallest  entries  were  correct  to  6  decimals.  Figure  1  shows  2_  and 
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z+  while  Figure  2  shows  the  bisectors.  An  unorthodox  representation  has 
been  used  in  order  to  emphasize  the  small  entries.  Instead  of  v{i)  we  plot 
(-1/  logio  |u(i)l)  sign(w(0)  and  the  vector  is  normalized  so  that  the  maximum 
value  is  1.  To  show  the  distortion  we  also  plot  (z+  +  z_)/\/2  in  the  standard 
manner  at  the  top  of  Figure  2. 

A  careful  look  at  z_  and  z+  shows  that  Wilkinson’s  idea  of  using  dis¬ 
joint  submatrices  could  not  be  satisfactory:  entries  z+(12), . . . ,  z^.(15)  and 
z_(5), . . . ,  z-(8)  are  not  negligible.  However  {z_,  z+}  is  the  most  sparse  ba¬ 
sis  of  the  invariant  subspace  for  Ai2)  A13.  So  there  is  no  basis  in  which  half  the 
entries  in  each  vector  are  negligible,  in  contrast  to  the  case  for  A205'^2i-  The 
overlap  of  the  submatrices  is  needed  to  obtain  the  middle  entries  accurately. 

These  remarks  do  not  contradict  the  fact  the  Wilkinson  constructed 
so  that  each  pair  of  eigenvectors  could  be  built  out  of  two  small  vectors  u 
and  V.  Here 


(1^1:10  _  ^  0,  -  Xl3l)v  =  0. 

Where  which  is  not  symmetric,  differs  from  by  replacing  entry 

(11, 10)  by  2.  The  eigenvectors  for  A12  and  A13  are 

(K(l),...,n(10),  0,  -u(10),...,-n(l))‘, 

(u(l),...,u(10),  u(ll),  u(10),...,u(l))‘. 

Wilkinson’s  idea  of  using  (1:11)  and  (12:21)  (with  =  0)  will  not  work  here 
because  u(7  :  10)  ^  v{7  :  10)  to  working  accuracy  in  contrast  to  the  case  of 

A205A21. 

Inverse  iteration,  even  with  well  chosen  right  hand  sides,  gives  poor  re¬ 
sults  unless  the  Gram-Schmidt  process  is  used  heavily.  The  output  from  the 
LAPACK  code  sstein  was  not  as  accurate  as  our  bisectors. 

The  following  sections  show  that  well  chosen  submatrices  yield  good  bases 
in  all  cases  when  the  cluster  is  reasonably  isolated  from  the  remaining  spec¬ 
trum. 
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6  The  Envelope  of  the  Invariant  Subspace 

Let  be  a  set  of  eigenvalues  of  T  well  separated  from 

the  rest  of  the  spectrum.  Let 

Tzi  =  Zi\i,  l|z,-l|2  =  l, 

and  define  I  :=  [A;,A/+m-i]  and 

Sj  =  span{zi, . . . 


The  envelope  vector  of  Sj  is  Sj  given  by 

€{j)  =  Sx{j)  :=  max{t;(i) :  u  G  5i,  llu||2  =  1}. 

The  extremal  vectors  in  Sj  are  characterized  by 

=  ||v''’ll2  =  l,  J  =  landn. 

We  consider  the  case  m  =  2  (close  pairs)  here  although  some  of  the  results 
may  be  extended  to  larger  clusters.  The  subvectors  that  are  useful  in  com¬ 
puting  an  orthogonal  basis  for  Sx  may  be  understood  as  approximations  to 
these  extremal  vectors.  A  little  more  notation  is  needed  before  the  results  of 
this  section  can  be  described.  By  Lemma  2  (proved  below) 


and  may  be  expressed  as 

=  zi  cos  ip  -f  Zi+i  sin  p, 


tan  p  = 


^/(l)  ’ 


and  p  plays  an  important  role  in  this  section.  We  may  assume  that  z/(l)  >  0, 
2/+i(l)  >  0.  The  Rayleigh  quotient  of  y^^^  is 

=  \i  cos^  p  +  A/+1  sin^  p.  (12) 


The  residual  of  y^^^  and  the  residual  norm  are 

3  y^^^  —  y^^ Vi  ~  Zi  cos(/5  (A/ —  pi) -b  2/+1  sin(^  (A/+1  —  pi), 


y^^Vi 

1/2  :=  IITy^i) -y(i)pi||2  =  (Ai -pi)^cosV  +  (A/+i -/>i)^sinV 

Ai+i  —  A/ 


=  (  sin2(p 


by  (12). 


(13) 
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defining  the  residual  norm  i^i  which  occurs  throughout  this  section. 

Our  goal  is  to  find  an  index  j,  1  <  j  <  n,  such  that  the  eigenvector 

e  J,  satisfies  =  0(A,+i-A,).  The  connection 

of  such  j  to  the  envelope  vector  €  is  given  by  the  principal  results  of  this 
section  which  we  summarize  first. 

For  each  j  such  that  has  an  eigenvalue  9^-^  €  I,  Theorem  1  shows 
that 


A  |s'‘'W)K0  + 1) « 


A<+i  ~  A/ 
2 


1/2 


« j/i, 


provided  that  (A/+i  —  Xi)fgap  1,  where 

gap  =  min{A/+2  —  —  A/_i}.  (14) 

For  and  the  dominant  pair  A20,  A21  one  finds  5(1)  =  5(21)  «  0.78  and 

5(2)  =  5(20)  «  0.58.  Thus  j  =  19  or  20  is  a  good  choice. 

If  does  not  have  an  eigenvalue  in  I  then,  by  Lemma  6,  "  has  an 

eigenvalue  in  I  and 

^,+1  |5^+^="(j  +  2)|5(j  +  l)« 


sin2^/’ 


A/+1  —  A; 


+  (0^+"="  -  Pnf 


where  =  zi  cosi/j  —  sin0,  pn  —  ■  Ty^ 

Note  that  ||5||2  =  and  the  average  value  of  5’s  entries  is  yf.  Theo¬ 
rem  1  tells  us  to  locate  5’s  maximal  entries  to  get  small  values  of  |•s^■•^(i)|. 

It  turns  out  that  the  vector  ^  ^ ,  for  suitable  j,  is  a  cheap  approximation 

to  y(^).  Throughout  this  section  we  abbrevaiate  x*’’’  by  Xj-  For  any  j  that 
keeps  \9^'^  -  /3i|/(A/+i  -  A/)  small,  Lemma  8  says 


-f  1)  = 


Xt^Pi)  .  1  2 


(  ( ~ 
\\  9<^P 


whereas 


s(^■)(^•-Kl)  =  I 


x.(0ib)^<'Ki)/(A---A), 

0, 


i  <  i, 
i  >  j- 
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Note  that  6^'^  only  occurs  in  the  second  (and  small)  term  in  +  1),  not 
the  dominant  first  term.  The  difference  between  and  depends  on  the 
two  ratios,  (A,+i  -  Xi)/gap  and  -  Pi|/(A/+i  -  A,).  In  the  same  vein  we 
show  (following  Lemma  7)  that,  for  the  appropriate  j  values. 


s(^)(l)  =  €{1) 


1  +  0 


( •^<+1  ~ 
V^ap(0i=-') 


if  5(1)  is  not  too  small, 


where  gap{6^''^)  =  min  \0  —  0^'^  over  all  eigenvalues  $  of  other  than 
and  is  slightly  different  from  gap  defined  in  (14). 

Indeed  very  good  approximations  to  the  eigenvectors  z/  and  Zi+i  are  given 
by  the  internal  and  the  external  bisectors  of  and  The  extremal 

vectors  and  are  not  quite  orthogonal: 


yW  .  y(")  «  (A,+1  -  A;)  ^  (A;  -  Aj)  ^ 


Proofs 

To  establish  all  these  results  in  a  simple  way  a  sequence  of  lemmata  will 
prove  useful.  Recall  that  Sx  =  span{zi, . . . ,  z/+m-i}- 

Lemma  2 


/-fm-l 

J(;)“  = 


1  =  / 


Proof: 


S{jY  =  maxu(j)^ 

1 

V  i=i 


=  max  < 
7 


=  1 


—  l  /+m  — 1 


/+m-l 

^  max  (  'yi  •  Y  j  =  S 


Equality  is  attained  when  ^  is  a  multiple  of  [zi{j)^ . . . ,  zi^rn-i{j))* 


Corollary  1 


n  /+m-l  l+m-1  n  (+m-l 

iifiP  =  E  E  E  E"‘0)'=  E 

jzzl  »=/  »=/  j=l 

Lemma  3  With  pi  =  •  Ty(^)  =  A,  cos^  9?  +  A/+i  sin^  y?,  then  for  any 

=  (A; -O^cos^¥’+ ('^'+1 

vi  =  sin2v7-l - 2 - ]• 

Proof. 

i/(^)  =  (A/  -  Pi  +  pi-  (f  cos^  (p  +  (A/+1  -  pi  +  pi-  0^  sinV 

=  j/j  +  2{pi  -  ^  +  -^'+1  “  /^i]  +  (^1  ~  0^ 

=  i'l  +  iPr-0^- 


Lemma  4  For  any  polynomial  x  of  degree  d  with  an  isolated  zero  6 


xiO  =  (e  -  0)x'{0) 


1  + 


1{±J) 

2  gap{9) 


(  9av{9)x"{9) 

\  X'{0) 


where  gap{9)  equals  to  the  separation  of  9  from  the  rest  of  x’s  zeros.  Also 


9ap{0)x"{0) 

X'(^) 


<  2{d-l). 


Proof.  A  Taylor  series  expansion  of  x  around  9  yields 


fn\  I 


1 


2.7/, 
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for  some  rj  in  the  convex  hull  of  (  and  6.  There  is  a  simple  expression 
for  x"{^)lx'{^)  if  denote  the  zeros  of  x  by  and  6  =  6j.  By 

logarithmic  differentiation, 


x'V,) 


=  2j2(e,-0tr' 


whence 


9api^j)x"{^j) 

xVj) 


Oj  -  0j±i 
Oj  —  6k 


<  2(d-  1). 


□ 


Remark.  For  evenly  spaced  zeros  the  upper  bound  is  2  ln(<i  —  1)  not  2[d  1). 


Lemma  5  For  each  j,  I  <  j  <  n,  either  ^  has  a  zero  in  [A/,  A/+i]  or 

Proof.  The  characteristic  polynomial  of  the  submatrix  of  T  obtained  by  delet¬ 
ing  row  and  column  j  is  By  Cauchy’s  interlace  theorem 

one  of  these  polynomials  (at  least)  has  a  zero  in  [A,', A,^|-i]  for  i  =  1. 

It  can  be  shown  that  if  (and  only  if)  both  polynomials  have  a  zero  in  [A,-,  A,+ij 
then  it  is  either  A,-  or  Ai+i. 

Now  we  can  establish  the  principal  result  of  this  section. 


Theorem  1  With  the  notation  developed  in  this  section,  if  Xj  has  a  zero  9^'^ 
in  (A/,A/.^i)  then 


+ 1) 


\9ap{e^-^)J  } 


Proof,  By  Lemma  2  and  (4)  in  Section  3 

S{j  +  iy  =  2^(i  +  l)'  +  ^^+l(i  +  l)^ 
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Square  Lemma  4,  for  x^''^  =  Xj  ^ 


1  +  0 


Az^Y 

_gap{e^-^) )  _ 

A/+1  —  0^' 


+  -  e'”f  [1  +  0 


and 


/?! 


I3,e{j  +  i)  =  \x'i(«'"’)\ 


1  +  0 


'A/+1  —  Af 

^gap{0^-:>)^ 


1^0 


By  Lemma  3, 

k 

To  relate  5(j  +  1)  to  the  submatrix  and  0^-^'s  eigenvector  recall  (6) 
in  Section  3 

s^i){l)s^^'>{j)x'ji0^-^)  =  l3f-/3j.i. 

Substitute  this  expression  into  the  equation  above  and  cancel  the  nonzero 
derivative  to  obtain 

I3j  ^(i  +  1)  =  ^  (^ap(0ib)) 


□ 


as  claimed. 

Lemma  6  For  all  j  such  that  €{j  +  1)  is  nearly  maximal  and  Xj  has  a  zero 
in  [A;,A/+i] 

=  €(iy  +  ■  rv  {«'■’)  sign  (^<^>(1))  1  +  0  ’ 

gap{j)  =  mm{Xi+2-0^'\0"-^  r  =  is  given  in  (15). 
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Proof.  To  analyze  in  terms  of  T  rather  than  note  that 
{T  -  ^ 

Use  r  =  ZAZ‘  to  obtain 

(  ““  ^  =  Z  (A z‘e,+ift»“’(i) 


and 


,»)(!)  .  ^ 

k=l 


is  taken  positive.  The  terms  corresponding  to  ^  =  /,  /  +  1  usually  dominate 
the  sum  and  it  may  be  written 

,«)(,)  =  +  ‘>  +  +  ^,,b,y)  (15) 


and 


=  E 


A:#/,(+l 


A/+1 

Zk{l)Zk{j  +  1) 
A,  -  ■ 


The  key  fact  is  that  the  first  two  terms  in  (15)  combine  to  give  5(1)E  Use 
(15)  to  obtain 

zi{\)zi{3  +  1)  ,  zi+i{l)zi+i{j  +  1)  _ 

X,  _  01-.:  A/+1  - 


=  «(1) 


2  x;(«'^o 


Pi'"  Pj 

Now  (6)  and  Theorem  1  give 

\Pi"'P, 


1  +  0 


( *^^+1  ~ 
\gap{e^-i) 


^(1) 

^(i  +  1)' 


(16) 


(17) 
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Multiply  (15)  by  and  then  use  (16)  to  find 


5(^)(1)2  =  £{l) 


,^/(^io>0)(l)3(i)(j) 


Now  use  (17)  to  simplify  both  terms 

A;+i  —  A/ 


1  .  (9  ^ 


sU){if  =  SiXf 


1  +  0 


gap{e^-i)j\ 


t[j  +  1) 


1  +  0 


\gap{e 


For  well  chosen  j  the  first  term  ^(1)^  usually  dominates  the  second.  When 
5(1)  is  small  then  j  must  be  large  to  ensure  that  6^'^  lies  in  [A/,Aj+i]  and  it 
is  possible  that  5(1)^  always  dominates  but  we  have  not  proved  that  yet. 
Note  that 

(i)  (^'=0  •  |t|  <  Efc#/./+i  \^k{i)zk{j  +  I)k IdMj)  <  ^  /9Mj) 

is  assumed  to  be  very  small. 

(ii)  £{j  +  1)  >  (above  average,  by  choice  of  j). 

(iii)  By  Lemma  7,  (i)  and  (ii),  whenever 

5(l)5(i  +  1)  >  2i/  /gap{j) 


then 


s(^)(l)  =  5(1) 


/  i/T  sign  {s^^'>{j))\ 

s{i)e(j  +  i)  J 


1  +  0 


(  ^<+1  ~  A/  \ 

\gap{6^--^) )  ^ 


i.  e.  s^-’^l)  «  ^(1)- 
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In  Figures  3  and  4  we  show  typical  envelopes. 

Figure  3  shows  snapshots  of  the  envelope  of  a  cluster  of  108  close  eigen- 
values,  each  differing  from  its  neighbors  by  0(c  •  spread),  from  a  matrix  of 
order  n  =  2053  supplied  by  George  Fann,  Pacific  Northwest  Laboratory, 
Richland,  WA.  The  matrix  arose  in  the  application  of  the  self-consistent  field 
Hartree-Fock  method  for  solving  the  nonlinear  Schroedinger  equation  for  Ze¬ 
olite  ZSM-5.  The  top  shows  the  first  50  entries  in  €  and  the  bottom  shows 
entries  470-550.  The  cluster  is  completely  determined  by  submatrix  (1:515). 
The  supports  of  all  the  eigenvectors  belonging  to  this  cluster  lie  in  (1:515). 

Figure  4  shows  the  envelope  of  span[zi2, 2^13)  for  Section  5.  The 

eigenvalues  are  close  to  6. 

The  ‘humps’  in  most  envelopes  are  unimodal  but  this  case  is  an  exception. 
As  shown  in  Section  5  either  peak  (4  or  6,  16  or  18)  will  serve  for  choosing  a 
submatrix.  However  it  is  vital  to  realize  that  indices  4  and  6  belong  to  the 
same  hump. 


Extremal  Vectors  in  Sj 


Lemma  7  For  any  0  €  (A/,  A;+i) 

(A  ■  ■  ■  “  Xi(/>i)  +  +  O  ((A|+,  -  A|)^)  , 

=  y"^-Ty>'K 

.  ci  A/+1  A/ 

Ui  =  sm  2(p  - - - - ■ 

Proof.  From  the  beginning  of  this  section 

=  zi  COS  zi+i  sin  tp,  (f  e  [0,7r/2], 

By  (4)  in  Section  3 

+  1)  =  .  cos  ^  •  Sin^. 


•  •  •  A 


Since  cost^  =  zi{l)/S{l)  and  5(1)  —  y^^\\), 

j/(i)(f  +  1)  ^  x.(A/)cos^y’  +  X«(A;+i)sin> 
y(i)(l)  “ 
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Expand  x.'  about  any  6  £  {Xi,  A/+i)  and  use  Lemma  3  to  find 


Xt(A/)  cos^  +  x«('^/+i) ~  (^  +  sin  (p) 

+  ‘f’  “  ^)  +  'f  {^i+i  “  ^)] 

+  ^x:'W‘'W'+o((Aw-A,f). 


The  coefBcient  of  xK*)  j“s‘  ft  “  ^  (P‘  “ 

right  side  simplifies  to 

Xi(pi)  +  ^  ((^>+^  ~  • 

as  claimed.  ^ 


For  comparison 


+  1) 
s(i)(l) 


f  Xi{0^'-^),  i  <  h 
\  0,  i>  j. 


The  closer  is  to  pi  the  closer  is  ^  ^  to  Let  us  consider  the 

magnitude  of  +  1)  which  we  approximate  by  0.  By  (6)  in  Section  3 


_  1 


(18) 


Hence,  with  6  =  9^'^  in  Lemma  7 , 

8/<‘l(j  +  1)  ^  ft  - 

y(0(l)  sP^(l)sPI(j)/?j 

1  I/,  ftx;(«’'')  sapjo'-’h'ne''’) 

+  2jnp(«W)'  ■  x'iif'J) 

+  higher  order  terms. 

By  Theorem  1,  for  the  best  values  of  j, 

0jS^^\j)£{j  +  1)  «  {9^'^)  =  [i^l  +  (^1  -  9^'^)'^]  ^  . 
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Recall  that  y(^)(l)  =  5(1)  and  5(1)  <  5(^)(1)  (but  «  if  5(1)  is  not  small). 
Thus 


y(l)(j  +  1)  = 


1 

+  t: 


Vi 


e(i  +  1) 


5(1) 

+  higher  order  terms] 


<  5(i  +  l) 


where 


M2  = 


[  V  (0i=i)  2  gap{d^-i) 

<  2(i  - 1) 


+ 


(19) 


pi  -  ,  ^1 

V  (0^=-')  gap{9^-^) 


gap{e^--^)x'l{e^-^) 
xW^) 

The  relation  (19)  tells  us  that  it  is  essential  that  both  ratios 

(Af+i  -  A/)sin2v? 

2gap{Q'^'j) 

be  small. 

Next  consider  Recall  that 

=  zi  cos  <p  +  2/+1  sin  V?,  =  z/  cos  V’  -  Zi+i  sin  V’, 

where  €  (0,7r/2).  This  choice  of  sign  for  y^"*  is  deliberate  and  yields 
€  (0,7r/2).  Since  z/  •  z/+i  =  0 

y(l)  .  y W  =  cos  cos  (p  —  smip  sin  <p  =  cos(V’  +  <p)- 

By  (6)  in  Section  3, 

z/(l)^/(«)x'(A/)  =  /?!  •  ■  •  0n-i  =  z;+i(l)2m(«)x'(A/+i) 

where  x  =  X^'"-  Hence 

/  _  -■^f+i(?^)  _  -^/(l)  X'(A;) 

zi{n)  z;+i(l)  x'(A/+i) 


1  (A;  —  A/+i)  ~  '^») 

tan  <p  (A/+1  -  A/)  nr(-^'+i  “  ^*) 
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where  []"  =  n"=i, </;,)+! 


.  So 


tan  (/?  •  tantp 


1  —  tan  (f  •  tanxj^ 
gap 


tt 

1 + (Aw  -  A,)  x;  ^ 


min|A,-  —  A;|,  i  ^  IJ  +  1. 


This  establishes 

Lemma  8 

y(i) .  y(")  =  COS  Ip  COS  (/>  (1  —  tan  (p  •  tanxp)  =  0 


7  The  Submatrix  Theorems 

We  consider  the  case  when  a  pair  of  adjacent  eigenvalues,  call  them  A_  and 
A+  are  much  closer  together  than  they  are  to  any  othei  eigenvalues. 

We  claim  that  there  is  a  submatrix  T^'-^  with  two  properties 


(a)  It  has  a  well  isolated  eigenvalue  6  in  (A_,A+). 

(b)  The  normalized  eigenvector  s  for  6  has  the  property  that  ^  ^  ^  is  very 

close  to  the  invariant  subspace  (under  T)  associated  with  A_  and  A+. 


If  we  apply  this  property  to  the  trailing  submatrices  Tfc,„,  for  some  k,  we 

^  ^  that  is  also  very  close  to  the  desired  invariant 


t 


obtain  another  vector 

subspace  for  A±.  Thus  ^  ^  ^  and  J  ^  make  an  excellent  basis  for  this 
subspace. 

The  smaller  is  A+  -  A_,  the  smaller  is  the  dot  product 
We  now  prepare  to  make  precise  the  proceeding  claims. 


f  s\  /  0  \ 

V  0  J  ■  V  t  / 
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Paige’s  Persistence  theorem  (in  Section  3)  says  that  there  is  at  least  one 
Ritz  value  in  the  Ritz  interval 


jl'>  ;=  |9!'>  -  ftWi,  +  ftwi),  Ui  = 


(20) 


for  all  /  =  1,2, ...  ,n  -  j.  When  /Jjo;,-  is  relatively  large,  as  happens  when  j 
is  small,  then  may  contain  several  eigenvalues.  However  as  a;,-  decreases 
then  becomes  disjoint  from  all  other  and,  as  j  increases,  shrinks  onto 
a  single  eigenvalue. 

Our  interest  now  centers  on  those  rare  cases  when  X-^  is  isolated  from  its 
neighbors  and  small  and  yet  2  eigenvalues  get  into  Xl^K  How  close  can  they 
be  to  each  other?  ^Ve  will  show  that  acts  as  a  barrier  in  the  sense  that 
the  two  values,  A_  and  A+,  satisfy 


A+  -  A_  >  Oi/3,uJi). 

f  5O)  \ 

So,  if  A+  -  A_  is  tiny  then  the  vector  I  q  j 
since 

\\T 


has  a  small  residual  norm 


Ji) 


0 


*1)'’ )  «PII  = 


More  notation  is  needed  to  state  Theorem  2. 


Secular  Equation 
For  some  j  <  n,  write 

/  ri:i  ,  \ 

•  Oj+i  •  ,  •  =  or  ^j+i. 

\  / 

If  (  is  not  an  eigenvalue  of  nor  of  then 

det{T  -  (1)  =  det{T^'-^  -  CI)[aj+i 

-  Cir^e,] .  -  O]- 

Consequently  the  eigenvalues  of  T  that  are  not  eigenvalues  of  T'^'^  nor  of 
j'j+2--ri  satisfy  the  nonlinear  secular  equation 

<7(A)  :=  a,+i  -  A  -  -  A/j'^e,-  -  -  A/)-'ei  =  0. 
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The  middle  term  will  be  written  as 


B^iof 


+ 4’\o 


where  Wi  =  =  |5{^^(j)l  and  €  (-^-i  A+),  and 


(21) 


and  the  final  term  as 

l3]+ie\{Tj+2,n  -  C^)"^ei  =  Tj+2,„(C)-  (22) 

The  following  results  bound  the  residual  norm  in  terms  of  the  sepa¬ 

ration  of  from  the  other  zeros  of  Xi-  Sometimes,  these  gaps  are  greater 
than  separation  of  from  eigenvalues  A  other  than  A±. 

Theorem  2  (Double  Occupancy)  Let  r„  denote  a  symmetric  unreduced 
tridiagonal  matrix.  Let  A_  and  A+  be  two  adjacent  eigenvalues.  The  notation 
developed  in  this  section  is  in  force.  Consider  those  indices  j ,  I  <  j  < 
that  satisfy  the  hypothesis 

(H)  There  is  a  single  Ritz  value  dy’  in  the  open  interval  (A_,  A+)  and 
neither  end  point  is  a  Ritz  value. 

For  such  j  test  the  condition  for  double  occupancy  of  in  terms  of  (21) 
and  (22): 

(VO)  V’P(^-)  +  'rj+2.n(A_)  <  Oj+i  -  <  V’/  (A+)  +  Tj+2,n(A+)- 

If  (VO)  does  not  hold  (only  one  eigenvalue  in  I^^^)  then 
-  X- > 

If  (VO )  holds  (both  eigenvalues  in  l[^^)  then 
>  A+  -  A_  >  +  2G, 


where 

G  =  G,(A_,A+)  =  i 


V’. 


^^^(A+)  -  V^P(A-)  ,  Tj+2,n{^+)  ~  'r?-t-2.n(^-) 


A+  -  A_ 


+ 


A+  -  A_ 
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Remark  1  Theorem  2  does  not  give  lower  bounds  on  A+  -  A_  because  G  is 
a  function  of  A_  and  A+. 

Remark  2  Gj  =  ^  + ’i+2,n(»?2)}  for  some  values  t/i  and  r]2  in 

(A_,A+).  Moreover  and  Tj^2,n  s^re  rational  functions  that  are  mono¬ 
tonic  increasing  between  poles.  By  (H)  (A_,  A+)  is  between  poles  of  ^.nd 
also  (A_,A+)  is  between  poles  of  Tj+2,n-  Hence  Gj  >  0  on  (A_,A+).  Our 
interest  is  in  those  j  for  which  Gj  is  smallest. 

Proof  {oi  Theorem  2).  Partition  T  -  (I  into  a  3  x  3  structure  with  blocks 
2^  .  j  j  1  ;  j  -|-  i^j  -f-2  :  n  and  then  form  the  Schur  complement  cr(C)  of  entry 
(i  +  l,j  +  !)•  Thus  r  -  C-f  is  congruent  to 

T^'-^  -(!&  a{C)  0  -  (I 


where 

<7{0  =  aj^^  - 

<T  =  0  is  sometimes  called  the  secular  equation. 

Observe  that  Ti=i0r^'+2:n  jg  submatrix  of  T  obtained  by  deleting 
row  and  column  j  +  1.  By  Cauchy’s  Interlace  theorem  the  Ritz  values  from 
and  interlace  the  eigenvalues  of  T.  By  (H)  (A_,  A+)  contains  Bf 

and  so  cannot  contain  any  Ritz  value  of  By  Theorem  5  in  [5]  A±  can 

only  be  Ritz  values  if  both  and  a  Ritz  value  of  coincide  at  A±. 

Thus  (H)  rules  out  this  possibility  and  so 

[A_,  A+],  the  closed  interval,  contains  no  Ritz  values  of  ^23) 

The  set  of  indices  j  that  satisfy  (H)  includes  j  =  n  - 1,  by  Cauchy’s  Interlace 
theorem  and  so  is  nonempty.  When  j  is  too  small  then  it  is  a  Ritz  value  of 
J’j-i-2:n  lies  in  (A_,A+)  and  not  B\^\ 

Let  !</<;■  denote  the  Ritz  value  in  (A_,  A+)  and  separate  it  from 
the  other  j-level  Ritz  values.  Thus 

-  ar'e,  =  +  tf  '(C) 
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where  ^P(C)  is  defined  by  (21). 

Recall  that  (5^(1), . .  •  is  a  normalized  Ritz  vector  for  0^  and 

^  J})  =  \s\^\j)\.  The  final  term  in  <t(C)  is  T,+2,n(C)  and  is  defined  in  (22). 
By  (H)  A-  and  A+  are  not  eigenvalues  of  Tj.  By  (23)  they  are  not  eigen¬ 
values  of  Tj+2,n-  Consequently  A_  and  A+  must  satisfy  the  secular  equation 

a2,  ,(i)^ 

<7(A±)  :=  ttj+i  -  A±  -  -  - - i’?  (-^i)  -  7'i+2,n(A±)  =  0. 

To  simplify  the  analysis  that  follows  write 

0  =  0\^\  (jj  =  P  =  ^j,  01  =  aj+i,  V’i  =  t  =  Tj+2,n- 

Next  rewrite  the  above  equations  in  terms  of  positive  terms  \+-9  and  6-X- 
to  find 

(A+  -  0f  +  2E{\+  -0)  =  ^V,  (24) 

(0- A_)2  +  2F(0-A_)  =  /3V,  (25) 

where 

E  =  E{\+)  =  [i/’/(^+)  +  ^(^+)  +  ^  - 
F  =  F{X-)  =  [a-0-  ^;(A_)  -  t(A_)]/2, 

and 

E  +  F  =  [t^/(A+)  -  il^i{X-)  +  t(A+)  -  t(A_)]/2. 

By  Remark  2  both  V*/  and  t  are  monotone  increasing  in  (A_,A+)  and  so 

£;-t-F>0. 

The  quadratic  a;*  +  2Ex  -  =  0  has  two  real  roots  whose  product  is 

First  we  establish  the  case  when  {VO)  does  not  hold. 

If  £"  <  0  then  the  positive  root  is  the  larger  (in  magnitude)  and  must 

exceed  I3u!.  Hence 

X+-0>  I3u},  ^  -  A_  >  0 

together  imply 

A+  -  A_  >  I3ijj. 
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Similarly,  if  F  <  0,  then  9  -  >  l3u,  X+ -  6  >  0,  and,  again 

A+  -  A_  >  I3uj, 

as  claimed.  Differences  smaller  than  (Suj  can  occur  only  when  both  E  = 
F(A+)  >  0  and  F  =  F(A_)  >  0;  that  is  Condition  (VO)  holds.  This  is  the 
case  considered  in  the  remaining  analysis. 

From  (24)  and  (25) 

\+-e  =  ^V/{  +  F2  +  F},  (26) 

e-X.=  +  f2  +  fy 


Thus 


A_ 


1 


+ 


1 


/32a;2  +  E2  +  e  +  F 

To  simplify  this  equation  note  that  the  function 


(27) 


/ -  ,  \/x^  +  —  X 

f{x)  :=  {x  +  \/a;2  +  /?2a;2)  = - 

is  monotone  decreasing  and  concave  upward  (/”  >  0)  for  a:  >  0.  By  concavity 
and  the  positivity  of  E  and  F, 


m + f{F) 

2 


E  +  F 


in  other  words,  from  (27), 

A+  -  A_  ^  /(F)  +  /(F) 

2/?2a;2  2 

To  simplify  this  relation  write 


>/ 


E  +  F 


(28) 


where  G  =  G(A_,  A.^)  is  given  in  the  statement  of  the  theorem. 
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Let  r  =  and  multiply  (28)  by  /5w  to  find 

j3u 

,Jl3W  +  +  (^)  G 

_ 1 _ 

^l  +  {rGy  +  rG' 

=  s/r+i^-rG. 

Add  rG  to  each  side,  square  and  subtract  r^G^  to  find 

r"(l  +  2G)  >  1 

as  claimed.  ^ 

Remark  3  In  the  cases  of  interest  to  us  j^io  is  small  and 

G  K  i  me)  +  Ae))  ■ 

The  next  task  is  to  turn  the  inequalities  of  Theorem  2  into  lower  bounds  on 
Ly  obtaining  upper  bounds  on  G.  These  bounds  depend  on  the  gap 

between  0\^'^  and  the  other  Ritz  values  of  and  T-'+2=”. 

Theorem  3  Let  A_  and  A4.  be  two  consecutive  eigenvalues  of  a  symmetric 
unreduced  tridiagonal  matrix  r„.  The  notation  developed  in  this  section  is  in 
force.  For  each  j,  1  <  i  <  n,  that  satisfies  the  two  hypotheses 

(H)  There  is  a  single  Ritz  value  in  the  open  interval  (A_,  A+) 
(GAP)  2yjUJi  <  gap{l),  where 

gapil)  =  -  €'v  l«!'’  -  •  =  1,"  -  J  -  l}. 

therij 

either  (A_,  A-|_)  ^  :=  ^ 

in  which  case  A-|.  —  A-  >  ^jtoi 

or{X-,\+)  _ 

in  which  case  ,  A^.  —  A_  >  2l3jLOtl  y/l  +  2rj, 
where 

Tj  =  {y]{l  -  (^f)  +  ^|+i)/2[5ap(0  - 
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Remark  4  Tht  obstwations  at  the  end  of  Section  3  say  that,  on  average, 
Fj  is  0(1);  the  geometric  mean  of  the  0’s  is  less  than  the  geometric  mean  of 
the  gaps  \6i  —  6i\,i^l. 


Proof.  It  is  only  necessary  to  bound  the  quantity  G  from  Theorem  2  and  the 
same  notation  will  be  adopted. 

For  some  rj  €  (A_,A+), 


0W(At)  - 

A+  -  A_ 


-  V-, 

=  ^j{i  -  u,?)/w!''. 


(29) 


where  is  a  weighted  harmonic  mean  of  Since 

the  minimum  value  never  exceeds  any  mean 


>  inin(0p^  -  r?)2 
=  min{(0[|\ 


ri?}- 


There  are  two  cases.  If  (A_,  A+)  0  then  A+  —  A_  >  0ju>i  for  the  same 
reasons  given  in  Theorem  2.  In  the  contrary  case,  one  has 

V  ^  [A-,  A+]  C  X^. 


Now  (GAP)  guarantees  that 

min{^|^\  ^|i\)  >  min{0/+i  -  Oi,Oi  -  0;_i}  -  0u)  >  0u. 

Similarly,  for  some  C  €  (A_,  A+) 


F?'+2,n(A+)  —  Tj4.2,n(A-)  _ 


A+-A_ 


I  fl(j+2.n 


„0+2,n) 


(0 


-  cj 


(30) 


where  is  a  weighted  harmonic  mean  of  —  C)^} 

and  so  (GAP)  yields 

7^0+2, n)  >  min(^P'’”^  -  C)' 

t 

> 

The  definition  of  gapil)  is  made  so  that 

>  {gapil)  -  >  {gap{l)  -  pujf , 

and  thus,  using  (29)  and  (30), 

2G  = 

^1(1  -  u;f )  +  PU 

{gap{l)  -  ^jW/)2  ’ 

:=  2r,- 

and  application  of  this  inequality  into  the  second  inequality  of  Theorem  2 
establishes  Theorem  3. 

Remark  5  When  we  first  established  Theorems  2  and  3  we  were  concerned 
that  the  gap  quantities  involved  Ritz  values  9]'^  etc.  and  not  eigenvalues  A,. 
However  the  study  of  some  challenging  examples  has  shown  that,  in  fact,  the 
gaps  involving  Ritz  values  do  lead  to  an  important  extension.  In  examples 
with  100  or  200  eigenvalues  in  a  cluster,  each  separated  from  its  neighbor 
by  about  10  or  20  ulps  (units  in  the  last  place)  combine  to  form  a  cluster 
of  nonnegligible  width  (e.  g.  1500  ulps).  We  find  that  for  each  submatrix, 
taken  with  its  nearest  neighbor,  the  gaps  of  Theorem  3  are  vastly  greater 
than  quantities  such  as  —  A/_i  and  A/+2  —  Oj  .  In  other  words.  Theorem  3 
can  be  applied,  submatrix  by  submatrix,  to  clusters  whose  total  width  is 
greater  than  0{e  •  spread). 

Recall  that  the  envelope  bound  in  Section  6  shows  that  there  are  ^-values 
for  which  our  gap  quantities  Gj  satisfy 
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In  fact  each  hump  in  an  envelope  contributes  1  to  H^H^,  on  average,  and  if 
the  hump  is  spiky  S{j  +  1)  ^  1  while  if  four  consecutive  entries  contribute 
to  the  hump  then  the  greatest  of  the  four  exceeds  1/2. 

8  The  Overlap  Theorems 

We  need  expressions  for  the  entries  of  the  vectors  constructed  to  span  the 
invariant  subspace  corresponding  to  a  close  pair  of  eigenvalues  A_  and 
of  r  =  r^=”.  Recall  that  x'''”(C)  is  the  characteristic  polynomial  of  the 

submatrix 

The  index  m  (not  unique)  is  determined  so  that  has  a  well  isolated 
Ritz  value  0^'"'  in  (A_,A+)  whose  normalized  Ritz  vector  has  a  small 
last  entry.  More  precisely 

Trailing  submatrices  are  used  to  find  a  suitable  index  I  (not  unique)  so  that 
has  a  well  isolated  Ritz  value  0‘'-^  in  (A_,A+)  whose  normalized  Ritz 
vector  has  a  small  first  entry.  More  precisely 

A-il/”(OI  =  0(A+-A_). 


The  spanning  vectors  are 


and  the  indices  1  <  /  <  n,  1  <  m  <  n,  play  a  crucial  role  in  the  analysis.  To 
avoid  index  troubles  the  entries  in  "  are  labelled  from  I  to  n,  not  from  1 
to  n  —  /  +  1. 

There  are  cases  in  which  m  <  I  and  then  p  and  q  are  orthogonal  because 
their  supports  are  disjoint.  So  we  turn  to  the  other  cases  when 

1  <  I  <  m  <  n. 

and  study  \p{j)qij)\  for  I  <  j  <  m.  We  will  show  that  p  and  g  are  nearly 
orthogonal.  There  are  formulae  for  the  magnitudes  of  p{i)  and  q{i) 

p{i)^  =  Xl,«-l(^)Xi+l.m(^)/xi,m(^)’  *  — 
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where  $  =  is  p’s  Ritz  value.  Similarly 

where  (p  =  is  q's  Ritz  value. 

It  is  readily  verified  that  for  tridiagonal  matrices 

Xl(i>n(C)  •=  Xl,t-l(C)Xt+l,n(C) 

is  the  characteristic  polynomial  of  the  submatrix  of  obtained  by  deleting 
row  and  column  i.  Hence,  by  Cauchy’s  Interlace  theorem  the  Ritz  values  of 
and  together  interlace  (weakly)  the  Ritz  values  of  In 

particular,  for  each  i,  the  open  interval  (A_,  A^.)  contains  either  a  Ritz  value 
of  or  a  Ritz  value  of  but  not  both.  In  an  exceptional  case  the 

open  interval  is  free  of  Ritz  values  just  when  Ritz  values  of  and 

coincide  at  either  A_  or  A+  (but  not  both)  and  in  that  case  p{i)q{i)  =  0. 

We  can  simplify  some  expressions  considerably  by  the  following  conven¬ 
tion: 

If  a  Ritz  value  G  [A_ ,  A+]  then  the  index  i  is  omitted.  Further 
denotes  the  smallest  Ritz  value  of  exceeding  A+  and  OL  denotes  the 
largest  Ritz  value  of  less  than  A_ . 


We  need  a  result  from  [3]. 

Theorem  4  With  the  notation  developed  above  let  Tsi  =  SiXi,  s*s,  =  1 
be  an  eigenvector  equation.  Let  denote  the  submatrix  obtained  from  T 
by  deleting  row  and  column  j;  its  spectrum  is  <  ...  <  Then  for 

I  <i  <n, 

.A.-e  ^p-A. 

'  ^  Aj  —  A,_i  Ai+i  —  Xi  Aj+2  —  Xi 
The  bound  for  i  =  1  and  i  =  n  is  obtained  by  omitting  quotients  with  out-of- 
bounds  indices. 

Proof.  For  tridiagonal  matrices,  see  [4], 

....  xfa>(Ao 

'  X'(A,) 
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By  Cauchy’s  Interlace  theorem 

Afc-l  <  ^  A/:  < 


So,  for  each  k  <  i,  the  quotient  <  1  and,  for  each  k  >  i,  the  quotient 

<1.  As  |ib  -  i|  increases  the  quotients  become  close  to  1.  By 
discarding  all  but  the  smallest  two  or  three  quotients  the  upper  bound  is 
obtained. 

In  what  follows  we  shall  apply  Theorem  4  to  submatrices  such  as 
and  One  of  the  principle  concerns  in  the  theorem  proved  below  is  the 

location  of  Ritz  values  such  as  or  9^^^^  .  By  (31)  9^^^  is  either  9^ 

or  and  we  shall  be  concerned  with  both  cases.  Recall  that  only  one 

of  and  lies  in  (A+,A++)  and  the  other  exceeds  A++,  the  next 

eigenvalue  of  T  greater  than  A+. 

Figure  5  is  worth  contemplating  before  reading  the  Overlap  theorem.  It 
shows  the  intervals  in  which  9\^~^  and  9^^^’^  will  lie  as  j  varies. 


Theorem  5  (Overlap)  Let  T  be  n  x  n,  symmetric,  unreduced,  and  tridi¬ 
agonal.  Suppose  that  adjacent  eigenvalues  A_  and  A^.  of  T  are  well  enough 
separated  from  the  remaining  spectrum  to  yield  the  indices  I  and  m  and  the 
vectors  p  and  q  described  at  the  beginning  of  the  section. 

For  each  j ,  I  <  j  <  m, 


where 

=  min{/5j_i,^j}, 

gap  =  Tain{gap{l),gap{m)}, 
gap{l)  =  mm{9‘-^-9‘l\9‘ff-9‘-'^}, 
gap{m)  =  min{0'=’"  - 

Proof.  First  confine  attention  to  those  j  values  such  that  G  (A_,A+). 

Consider  the  expression  for  p{j)^  in  Theorem  4  using  ™  and  extract  the 
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three  smallest  terms  in  the  product  to  find 


pOT  ^  \ 


I  0l-m^0l-rn  *  0Um^Ql'.m  9 


(32) 


otherwise. 


Thus  the  middle  term  itself  is  bounded  by  (A+  -  X-)l9ap{m).  A  smaller 
bound  emerges  by  considering  either  a  neighboring  term  from  p(i)^  or  the 
smallest  term  in  q{jy.  Without  loss  of  generality  we  suppose  that  the  closest 
Ritz  values  outside  (A_,  A+)  are  on  the  right.  By  Cauchy’s  Interlace  theorem 
the  open  interval  (A+,A++)  contains  either  0+^  ^  or  "  but  not  both. 

Case  1:  €  (A+,A++).  ,  /,i  .i 

To  obtain  a  bound  better  than  A+  -  A_  on  1^^™  -  1  consider  0  as 

an  approximation  to  and  apply  the  Gap  theorem  for  Rayleigh  quotients, 

r  ii-'p  I 


where 


gap{m, 

||r||  =  IKT'’™  -  e^-^-^lmYsW  = 

, = ( = \s'->-'u  -  i)i. 


(33) 


and 


gap{m,j  -  1)  =  min  {9]^^  -  -  ^1=”*} 

For  future  application  note  that 


gap{mj  -  1)  =  9ap{m) 


1  +  0 


V  gap{m) 


(34) 


since  and  lie  in  (A_,A+).  To  bound  ||r|l  we  apply  the  Double 

Occupancy  Theorem,  proved  above,  not  to  but  to  T.  Single  occupancy 
guarantees  that  ||r||  <  A4.  —  A_  and  hence 


,  (A+-A_2 

ml:m  _  0ly-l  <  ^  ± 

'  5ap(m) 


^A+-A_' 
1  +  G  (  r 

V  9m^) . 
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so  that  |p(j)|  <  (A+  -  \-)l gap.  This  is  already  tighter  than  the  bound  to  be 
established.  Double  occupancy  yields  a  weaker  bound: 

||rl|2<(H-2G,_i)(A+-A_)2  (35) 


where 


/i  =  -(A+  +  A_). 


1  +  0 


V  . 


Moreover,  from  Remark  3  in  Section  7  and  (29),  (30), 


xl>\p)  < 

t\p)  < 


(36) 

(37) 


In  Case  1,  >  t'  and  putting  (35)  into  (33)  yields 

A..-A_ 


^ - 7—\ — ^  < - 7~\ 

gap{m)  gap{m) 


(1  +  20j_i ) 


A+-A- 

gap{m) 


} 


Now  we  must  use  the  right  hand  term  in  (32),  noting  that 


0^)m 

~e^ 


- 

- 


^W-1  _  0l:m 
QUm  _  Ql-.m 


( -  p 

I  _  01  :m 


1+0 


A.f-A-\ 

gap{m)  ) 


■  (39) 


In  the  analysis  to  follow  we  shall  drop  the  1  from  1  +  20j_i  because  it 
contributes  only  a  higher  order  term  to  the  bounds.  Insert  (36)  into  (38)  and 
multiply  by  (39)  to  find 


P{3?  < 
< 


0\.m  _  0^”^  _ 


(A+-A_)  .  /  2?U  (A-,-A-)l 

jop(m)  "'”1  WW  J 
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oT'  -  M 

gap{m) 

(A+-A-) 

gap{m) 

i  +  o(— 


1  +  0 


p-h  -  A- 
V  5«P(”^) 


+  0 


p-^-A- 

V  5ap("i) 


mm 


A+  -  A-  I 

I  5ap(m)  ’  -  p  5ap(m)2  j 


-X. 


+  0 


pap(m) 

Bound  the  min  by  the  geometric  mean  to  find 


/A.^  -  A- 

V  gapi'^) 


(40) 


p{jf  <  -41 


A+  -  A-  ft-i 


'A+  -A_' 


1/2 


gap{m)  gap{m)  \gap(m) 

+  0 


^  (K  -  A- 

V  gap{m) 


\  gap{m) 


(41) 


Since  q{jY  <  1  the  claimed  bound  holds  in  Case  1. 


Case  2:  €  (A+,  A++).  _ 

The  argument  is  similar  to  Case  1  but  now  it  is  q{j)  that  offsets  a  large 

value  for  Gj-i.  In  Case  2,  r'  >  rp'  and  so  (38)  yields 


2D] 


A+  -  A_ 


1  +  0 


A_ 


i- 

\  gap{i'n) 


(^^+i:n  _  gap[m) 
+  0 


f  A+ 


I  pap(m)  , 


The  smallest  term  in  q{jY  gives 

ghn  _  Ql-.n 


gur 


-  p 
gap{l) 


1  +  0 


A+  -  A- 
gap{l) 


(42) 


Now  take  the  product  of  (40)  and  (42)  and  bound  min  by  the  geometric  mean 
to  find 


A+-A_ 


-  p  213] 


A+  -  A-  \ 


pij)\U?  <  gap{m)  [  gap{l)  ’  -  p  '  gap{m)gap{l)  j 
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< 


1+0 
A 


A+-A_ 


^/2 


gap 

4  -  A- 


A_ 


gap 


.rit. 


gap{m)  gap{l)  \  gap{m) 


A+-A- 

gap 


(43) 


We  see  that  (43)  is  an  instance  of  the  bound  claimed  in  the  theorem. 

The  analysis  for  j-values  in  which  €  (A_,  A+)  is  the  dual  of  what  is 
given  above.  Attention  concentrates  on  q{jy  instead  of  p{j)^,  go^pi})  replaces 
gap{m)  and  the  roles  of  and  are  exchanged.  By  Cauchy’s  Interlace 
theorem  either  €  (A_,A+)  or  6  (A_,  A+)  for  each  j,  1  <  i  <  n, 

and  so  the  proof  is  complete. 


Remark.  The  analysis  shows  that  the  configuration  of  Ritz  values  has  to  be 
quite  special  to  yield  a  value  for  lp(i)9(i)l  as  large  as  0((A+  -  \-)lgapfl'^. 
If  in  Case  1,  or  in  Case  2,  is  close  to  A+  or  to  then  the 

minima  in  (40)  and  (41)  are  0((A+  -  X-)lgap)  and  that  same  bound  holds 
for  |p(i)9(i)|-  Since  there  is  usually  a  little  freedom  in  the  choice  of  I  and  m 
we  can  expect  that  for  j  close  to  /  and  to  m 

IpOXj)!  =  0  ■ 


It  follows  that 


Overlap{p,q)  :=  \p\  -^1  =  0 


9  Accuracy  of  Subspaces 

Consider  the  basis  {p,q}  produced  by  the  use  of  submatrices.  Here  p’s 
support  is  1  :  m,  and  g’s  support  is  /  :  n  and  f,  m  are  chosen  so  that 
\l3mP{m)\  and  are  small.  In  addition  \\p\\  =  ||g||  =  1.  We  consider 

how  well  span{p,  q)  approximates  the  invariant  subspace  for  A_  and  A+. 

Let  =  p*Tp  be  p’s  Rayleigh  quotient  and  let  0'-"  =  q*Tq  be  q’s 
Rayleigh  quotient.  We  have 

Tp  =  pO^''^  +  em+i^mP{m),  Tq  =  q9‘''"- +  ei-i^t-iq{l).  (44) 
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There  is  no  loss  in  shifting  T  to  the  mean  /x  of  the  two  Ritz  values  0^  '^  and 
Also  let 


6  = 


0Um  _  0l:n 


There  are  two  expressions  for  k  =  p*{T  — 

K  =  p*q6  +  q{m+l)l3mp{m) 


(45) 


Also 


Then 


K  =  -p*q6  +  p{l  - 

b.«rb.«i  =  ( V)’ 

/  Ql'm  ^  \ 

\p,QYT\p,q]=  ^  0t:n  J- 

R:=  {T  -  pl)\p,q]-\p,q](^l  ) 


By  (44)  and  (45), 

R  =  [em+i0mp{'m)  -  qii,  ei-\Pi-\q{l)  -  P^]. 

Recall  the  supports  of  p  and  q  and  use  (45)  to  find 

0ip{^?  -  p-qK^  ^ 

^  ^  -  y  p*qK^,  -  2p*qSK  J 

Since  p  and  q  are  not  orthogonal  the  proper  measure  of  [p,g]’s  residual  is 


-1/2' 


Thus  <7^  is  the  largest  zero  of 


det 


RTR  - 


2/1  P’9 

p'q  1 
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So 

[1  _  {p*qf]<T*  -  +  Phqilf  -  2«^  -  2K\p^qf](r^  +  const  =  0. 

(Tmax  is  majorized  by  the  sum  of  the  roots 

,  ^  &M!±^wi(01z2£ 

^max  <  j  _  (P*g)^ 

Now  (45)  may  be  rewritten  as  2ac  =  £  +  A1,  defining  C  and  A4  in  a  natural 
way.  Since 

2^  =  C^  +  M^-\iC-Mf 

we  obtain 

«^max  <  1)^)  + ~  1)^) 

+  \[0mP{m)qim  +  1)  -  I3i-iq{l)p{l  -  1)]^}  •  [1  +  (p*9)^]-  (46) 

4U 

This  is  an  easily  computed  bound.  The  closer  is  \q{'m  +  1)|  to  ||9||oo 
\p{l  —  1)1  to  IIpIIoo  the  lower  is  the  bound  on  amax-  From  standard  gap 
theorems  in  [4]  the  sine  of  the  error  angle  is  less  than  cTj^ax/ 9^Pi  where  gop 
is  the  separation  of  [A_,A+]  from  the  rest  of  the  spectrum. 

The  General  Case 

Given  are  P11P21  •  •  •  iP#  with  the  support  of  Pj  on  {Ij  :  mj).  Also  p^  •  Pj  0 
if  li  -  i|  >  1  and  ||p_,|l  =  1,  j  =  1, . . . ,  #.  By  construction, 

Tpj  =  PjOj  +  Clj-\^lj-\Pj{lj)  +  Cmj+l^mjPji'fT^j)-)  Oj  =  0  ^  (47) 

In  order  to  simplify  expressions  it  is  convenient  to  shift  T  to  the  mean  value 
p  of  the  Ritz  values  ,  0#-  Let  0,-  =  /x  +  6,,  i  =  1, . . . , 

There  are  two  expressions  for  k,-  :=  p*{T  —  pI)Pi^i- 

_  f  Pi+i  •  Pi  +  Pt+l(^t ■  +  l)^miPi(^^i)5 

”  1  Pi  -Pi+l^t+l  +P.(/t+l  -  l)/5/,+i-lPi+l(Wl)5 

because  of  the  disjoint  supports. 

Define 

n  :=  {T  -  pI)Pi  -  SiPi  -  K,-iPi-i  -  ^iPi+i 


(49) 
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where  kq  =  k#  =  0.  Let  R  =  [pi, . . . ,r#].  By  (47), 

r,-  =  +  em,+i^m,Pi(^»)  -Pi+l^^i- 

Again  the  disjoint  supports  show  that 

Ti  •  r.-i-i  =  Pi_i  •  Pi  +  Pi  •  Pi+l  +  Pi+l  ’  Pi+2  '^i^i+1 )  (50) 

Ti  •  r,+2  =  K.«.+l 

and 

Vi  Vi  =  ^l.iPiihY  + 

-2Ki.iPi-x{li  -  l)/3(.-ip,(/,)  -  2/c,pi+i(m.  +  l)0m,Pi{^^i)- 

By  (48), 

K?  -  2KiPi+i{mi  +  l)l3miPi{^i) 

=  [/c,-  -  Pi+l{mi  +  l)^miPi{'mi)Y  -  Pi+l{^i  +  l)^/^miP«'(”^«)^ 

=  (Pi  •  Pi+l  ^iY  -  Pi+l  i^i  +  ^Y^lM^iY- 

Thus,  for  i  =  1, 2, . . .  #, 

Vi  •  Ti  =  (1  —  Pi-l{h  —  1)^)  0li-l  Pi{U) 

+  (1  -  p,+i(m.-  +  \Y)^lM'‘^iY 

+  {Pi-i  •  Pi  ^i-iY  +  {Pi  •  Pi+l  ^iY^ 

with  out  of  range  terms  set  to  zero  when  i  =  1  and  Let 

/  Pi -Pa  P2P3  •  P#-iP# 

/;=  tridiag  I  1  1  •  • 

\  P1P2  P2  P3  •  P#-i  •  P# 

Then  the  measure  of  the  quality  of  Range  [Pi,  ■  •  •  ,P#]  as  an  invariant  sub¬ 
space  is  given  by  CTmax  where  is  the  largest  zero  of 

det[R*R-al^/l].  (53) 

some  terms  in  (48),  (50),  and  (52)  are  much  smaller  than  others. 
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Thus  if  w  :=  cluster  width 

{R*R)ii  =  (1  -  Pi-i{h  -  1)")  AVi  Piihf 

+(1  - 

{R^Rki+i  =  0  +  O(u;''/') 

{R*R)i,i+2  =  KiK.+l- 

Using  this  approximation  it  is  straightforward  to  approximate  the  largest 
zero  of  (53)  by  bisection.  The  largest  diagonal  entry  of  R* R  is  a  reasonable 
approximation  to  <r^ax* 

10  Counting  Ritz  Values 

In  order  to  justify  the  selection  of  submatrices  some  background  material  is 
needed.  Recall  that  denotes  a  Ritz  value  from  T'bV,  the  submatrix 

obtained  by  deleting  row  and  column  j  from  T'"'". 

•  Cauchy’s  Interlace  theorem  (true  for  symmetric  matrices): 

Let  :=  For  each  j  =  1, . . . ,  n 

Ai  <  <  Ai„  < 

•  Unreduced  Tridiagonal  Interlace  theorem:  For  each  j  <  n,i  <  j,  and 
1  <  /  <  n  —  j  there  is  a  Ritz  value  in  the  closed  interval  [6]'\ 

For  a  proof  see  [3]. 

•  The  Window  Count. 

Let  I  be  any  fixed  closed  interval  on  the  real  line.  Let  #r(i  ;  k)  be  the 
number  of  Ritz  values  of  in  J.  Then  #j(i  :  is  ‘nearly’  monotone 

increasing  with  k.  More  precisely, 

-1  <  #i{j  -k  +  l)-  max  :  i)  <  1,  for  all  k. 

This  result  is  a  direct  corollary  of  the  Tridiagonal  Interlace  Theorem. 

Lemma  9  (Window  Count)  Let  1  he  the  convex  hull  of  a  cluster  of  ad¬ 
jacent  eigenvalues  of  unreduced,  symmetric,  tridiagonal  T  =  .  If  the 

window  count  ’■  j)  is  well  defined  for  _)  =  2, . . . ,  n  —  1  then 

#i(l  :;-!)  +  #r(i  +  1  :  n)  =  #r(l  :  7r)  -  1,  ;  =  2, . . .  ,n  +  1. 
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Proof.  Since  the  end  points  of  I  are  the  extreme  eigenvalues  of  the  cluster 
there  are  exactly  #t(1  :  n)-l  abutting  subintervals  in  J  of  the  form  [A,-,  Xi+i]. 
If  the  window  count  is  well  defined  there  are  no  zero  pivots  in  the  triangular 
factorization  (up  or  down)  with  shifts  at  the  end  points  of  I.  Hence  no  Ritz 
values  of  or  fall  at  the  end  points  of  I.  By  Cauchy’s  Interlace 

theorem  (one  can  assign  a  Ritz  value  of  to  each  subinterval)  there  are 

either  #i(l  :  n)  -  1  or  #r(l  :  n)  Ritz  values  of  in  J,  for  each  j.  By 

the  ‘coincidence’  property  of  tridiagonals  there  can  only  be  #r(l  =  ”) 
values  in  I  if  one  of  I’s  end  points  is  a  Ritz  value  and  this  is  ruled  out  by 
the  assumption  on  the  window  count.  ° 

11  Submatrix  Selection 

At  present  we  have  no  preferred  method  for  choosing  submatrices  automat¬ 
ically.  Below  we  present  two  methods  that  have  been  satisfactory  so  far. 


Mid-Point  Selection 

T  and  I  are  given.  I  is  the  convex  hull  of  a  cluster.  Let  #  =  #i(l  :  n). 
Define,  for  j  =  0, 1, . . . , 

Ij  =  max{z  :  #j(l  :  i)  =  j}-  (^4) 

Define,  for  j  =  1, . . .  ,#, 

rrij  =  + 

Note  that  7#  =  n.  Let  7_i  =  -1.  Take  as  initial  submatrices 

(ij-2  +  2  :  nij),  j  =  1, . . . ,  ^. 

Justification.  Since  #i(l  :  f)  is  nearly  monotone  increasing  in  i  the  in¬ 
dices  {fj}  are  strictly  monotone  increasing  in  j  thanks  to  the  max  in  their 
definition.  Hence 

^  7j  <  7j  -f  2,  j  =  0, 1, . . . ,  (55) 


Hence 


rrij-i  <  /j-2  <  h-'i  +  2- 
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Thus  the  supports  of  the  submatrices  are  disjoint  except,  possibly,  for  nearest 

neighbors.  •  c- 

Next  we  show  that  there  are  ‘enough’  Ritz  values  in  each  submatrix.  Since 

X  is  not  the  convex  hull  of  the  Ritz  values  of  the  submatrix  (fj-2  +  2  .  n)  in  X 
the  window  count  lemma  is  not  applicable.  By  Cauchy’s  Interlace  theorem 

#j(/,_2+2  ;  n)-l  <  #i(rj-2+2  :  mj)+#i(mj+2  :  n)  <  #r((,_2+2  ;  n)  +  l. 

(56) 

However  Lemma  9  may  be  invoked  twice  to  obtain 

^j(mj  +  2  :  n)  =  (^  —  1)  —  #j(l  ^  —  j  —  1,  (57) 

and  _ 

#j(/,_2  +  2  :  n)  =  (#  -  1)  -  #t(1  :  h-i)  =  #  -  i  +  1-  (58) 

Use  (57)  and  (58)  in  (56)  to  obtain 

3>#x(/;-2  +  2:m,)>l.  (59) 

In  the  exceptional  case  that  ^i(/j_2  +  2  :  mj)  >  1,  for  some  j,  the  support  of 
the  jth  submatrix  may  be  reduced  from  either  or  both  ends  until  the  count 

is  exactly  1.  .  ^ 

There  is  a  dual  algorithm  using  trailing  submatrices  that  delivers  mj, 

j  =  1  ,...,#+  1  such  that 

mj  =  min{i  :  #i{i  :  n)  =  #  -  j  +  1} 

and  mid-points  Ij,  j  =  1,...,#,  with  /,■  =  +  mi+i)/2j •  This  process 

yields  more  balanced  submatrices 

(Ij-.ruj),  ;  = 

with  the  same  bounds  and  disjoint  support  properties. 

We  used  a  selection  very  close  to  {(/j  :  mj)]  in  1989  before  our  analysis  of 
close  pairs  was  developed.  The  performance  was  very  satisfactory  but  there 
is  no  reason  why  the  mid-points  in  the  ranges  {Jj-i  +  1, . .  .,lj}  should  give 
the  smallest  coefficient  of  \I\  in  the  residual  norm  bounds.  We  want  the 
index  i  in  (/j—i  -j- 1, . . . ,  Ij}  that  gives  a  minimal,  or  small,  value  to  the  Gi  of 
the  Double  Occupancy  Theorem. 
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Selection  by  Pairs 

The  previous  results  for  close  pairs  of  eigenvalues  may  be  used  in  a  system¬ 
atic  way  to  produce  appropriate  submatrices  for  isolated  clusters  containing 
any  number  of  eigenvalues.  The  goal  here  is  to  show  the  existence  of  the 
submatrices,  not  to  produce  an  efficient  algorithm. 

From  Section  7  if  #(1  :  m)  =  2  then  there  exist  suitable  indices  //  and 
V  (by  no  means  unique)  such  that  #(1,^')  =  1  and  ^{n,m)  =  1.  Usually 
<  1/  but  that  is  not  necessary  in  what  follows.  There  are  basis  vectors  with 
supports  on  (l,i/)  and  on  whose  residual  norms  are  proportional  to 

the  separation  of  the  two  Ritz  values  that  cause  t^(1  •  w)  2. 

Now  suppose  that 

#;=  #(l,n)  =  #j(l,n)>2. 

Let  h  be  maximal  such  that  #(l,/i)  =  1.  For  the  submatrix  (1  :  h  +  1)  let 
the  optimal  indices  {fi  and  u)  be  72  and  ki.  The  Ritz  vector  (with  Ritz  value 
in  J)  for  submatrix  (1  :  fci)  is  the  first  basis  vector.  Set  ji  =  1.  The  Ritz 
vector  for  (j2  :  h  -|-  1)  is  not  used  but  the  index  j2  will  play  a  role.  Check 

that 


#(;2  :  n)  =  #(1  :  n)  -  1.  (60) 

If  not,  adjust  j2  until  (60)  holds.  Having  peeled  off  a  submatrix  ji  :  ki  from 
the  top  and  then  discarded  rows  1  :  j2  -  1  from  T  we  proceed  in  the  same 
way  at  the  bottom  of  T. 

Let  p  be  minimal  such  that  il^{p-,n)  —  1.  For  the  submatrix  (p  —  1  :  n) 
let  the  optimal  indices  be  j#  and  The  Ritz  vector  (with  Ritz  value  in 

I)  for  submatrix  (i#,n)  is  the  final  (#)  basis  vector.  Set  4  =  n.  The  Ritz 
vector  for  submatrix  (p  —  1  :  4-i)  %-i  P^ays  a  role.  Check 

that  , 

#(i2,4-i)  =  #  “ 

If  not,  adjust  4_i  until  (61)  holds.  If  #  -  2  >  2  then  repeat  the  proce¬ 
dure  just  described  on  {32  :  4-i)  to  obtain  two  new  submatrices  (;2,A;2), 
(i#-i5^#-i)  and’  possibly,  a  remaining  submatrix  with  fewer  Ritz  values  in 
I.  Eventually  one  obtains  #  submatrices  {ji  :  h),  i  =  1,...,#  each  of 
which,  by  construction,  has  a  simple  Ritz  value  in  J.  The  associated  Ritz 
vectors,  with  zeros  appended  to  make  n-vectors,  constitute  a  good  basis  for 
X’s  invariant  subspace. 
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Should  it  ever  occur  that  j,+i  <  ki-i  we  are  at  liberty  to  increase  ji+\  or 
decrease  a  little  subject  only  to  the  constraint  that  #(i,-  :  h)  =  1  for 
each  i.  Our  goal  is  to  have 

ki-\  <  max  entry  in  s-’’’’*''  <  j.+i,  i  =  2,  —  1. 

Selection  by  Envelope 

A  way  to  choose  submatrices  is  suggested  by  Section  6.  First  find  the  envelope 
vector  S  and  then  find  #  entries  of  S  that  are  local  maxima.  Suppose  first 
that  there  is  a  unique  set  of  such  positions  /ji,  ^2, . . . ,  k#-  For  1  <  j  <  #  the 
jth  submatrix  is 

(kj-i  +  1  :  kj+\  —  1).  (62) 

Usually  the  first  and  the  last  are  (1  :  /:2  —  1)  and  (^#-i  +  1  -  n).  However 
in  general  we  must  be  more  careful  and  define  ko  and  Let  the  first 

nonnegligible  entry  of  £  be  in  position  /jq  +  1  and  let  the  last  nonnegligible 
entry  be  in  position  ~  !•  Now  (62)  also  gives  the  submatrices  for  j  1 
and  j  =  #• 

In  case  the  location  of  the  zth  summit  is  given  by  several  indices  then  take 
ki  to  be  that  set  and  interpret  fc, •  +  1  as  1  +  max  ki  and  ki  —  1  as  —  1  +  min  ki. 
To  quantify  the  adjective  negligible  we  propose  a  threshold  of  macheps  •  ||5|| 
=  macheps  •  \/^. 

The  approximation  of  £  is  an  implementation  issue  that  will  not  be  ad¬ 
dressed  here. 


12  More  Examples 

A  Glued  Wilkinson  Matrix 

In  Section  5  we  studied  .  Here  we  use  but  take  4  copies  and  connect 
them  by  an  off-diagonal  entry  e  which  is  called  ‘the  glue’  in  the  matrix  VFioo- 
In  an  obvious  extension  of  our  notation  in  Section  2 

Wioo  =  tridiag  ^  VFjt  1^2'5  ^25  ^  • 
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If  e  is  too  small,  like  10-^  then  W^ioo  is  too  close  to  a  direct  sum  of  4  matrices 
and  calculation  of  orthogonal  eigenvectors  is  not  hard.  If  e  is  too  large  (>  2) 
then  the  eigenvalues  are  sufiiciently  well  separated  to  be  treated  as  isolated. 


and  show  the  largest  8  computed  eigenvalues  in  Table  2.  Without  the  use  of 
Gram-Schmidt  orthogonalization  inverse  iteration  gives  neither  small  resid¬ 
uals  nor  adequate  orthogonality. 


A93 

12.577864 

A94 

12.577870 

A95 

12.577881 

A96 

12.746191 

A97 

12.746193 

A98 

12.939114 

A99 

12.939115 

Aloo 

12.939117 

Table  2:  Selected  Eigenvalues  of  Wioo 


Figures  6,  7  and  8  show  the  submatrix  indices  used  and  the  basis  vectors 
they  yield.  The  vectors  are  plotted  on  a  logarithmic  scale  with  the  correct 
sign  attached.  All  entries  less  than  10“®  are  treated  as  0.  The  reason  for 
using  log  scale  is  the  to  focus  attention  on  the  smallei  entries. 

In  this  example  the  submatrices  overlap  by  only  one  or  two  indices.  The 
residual  norms  are  the  magnitudes  of  the  first  and  last  entries  and  are  all  less 
than  tncichcps  •  spread.  The  dot  products  between  vectors  in  each  group  are 
almost  zero  because  the  supports  are  almost  disjoint.  On  the  other  hand  the 
supports  of  the  vectors  X93  and  Xgg  are  identical  and  orthogonality  comes 
from  cancellation.  These  dot  products  are  less  than  30  •  vnacheps. 

An  Example  from  the  Lanezos  Algorithm 

The  Lanezos  Algorithm  with  no  reorthogonalization  was  run  in  double  pre- 
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cision  on  a  diagonal  matrix  of  order  205 

D  =  diag{l,2, . . . ,  200, 400, 400, 400, 400, 600) 
with  starting  vector 

The  run  stopped  at  step  87  and  the  resulting  tridiagonal  matrix  Ts?  = 
had  5  copies  of  600,  four  copies  of  400,  and  a  single  eigenvalue  at  200,  to 
single  precision. 

Figure  9  shows  the  four  vectors  corresponding  to  the  cluster  at  400.  All 
these  calculations  were  in  single  precision.  Note  that  the  overlap  of  the 
supports  is  greater  than  in  the  glued  Wilkinson  matrix.  However  all  nonzero 
products  were  about  10"^“*,  much  les  than  macheps. 
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Figure  7:  Basis  eigenvectors  of  VFioo  for  {^96)  -^97} 
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Figure  9:  Eigenvectors  0:79,  »80,  *81,  3^82  for  Tst  on  a  log  scale 
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